Given the EEA reference grid of Belgium at 1km resolution and the Belgian Natura2000 protected areas, this document assesses the cells covering the protected areas.
library(tidyverse) # To do datascience
library(tidylog) # To provide feedback on dplyr functions
library(here) # To find files
library(sf) # To work with geospatial data
library(BelgiumMaps.StatBel) # To load Belgian administrative boundaries
library(INBOtheme) # To apply INBO style to plots
library(leaflet) # To make interactive maps
We import UTM grid data of Belgium at 1 by 1 km resolution. All grids have CRS (Coordinate Reference System) equal to Belge Lambert 1972:
be_grid <- st_read(here::here("data", "external", "utm1_bel"))
## Reading layer `be_1km' from data source `C:\Users\damiano_oldoni\Documents\INBO\repositories\pipeline\data\external\utm1_bel' using driver `ESRI Shapefile'
## Simple feature collection with 51726 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 3768000 ymin: 2926000 xmax: 4080000 ymax: 3237000
## CRS: 3035
be_grid_crs <- st_crs(be_grid)
be_grid_crs
## Coordinate Reference System:
## User input: 3035
## wkt:
## PROJCS["ETRS89 / ETRS-LAEA",
## GEOGCS["ETRS89",
## DATUM["European_Terrestrial_Reference_System_1989",
## SPHEROID["GRS 1980",6378137,298.257222101,
## AUTHORITY["EPSG","7019"]],
## TOWGS84[0,0,0,0,0,0,0],
## AUTHORITY["EPSG","6258"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.01745329251994328,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4258"]],
## UNIT["metre",1,
## AUTHORITY["EPSG","9001"]],
## PROJECTION["Lambert_Azimuthal_Equal_Area"],
## PARAMETER["latitude_of_center",52],
## PARAMETER["longitude_of_center",10],
## PARAMETER["false_easting",4321000],
## PARAMETER["false_northing",3210000],
## AUTHORITY["EPSG","3035"],
## AXIS["X",EAST],
## AXIS["Y",NORTH]]
GIS data of protected areas as defined in Natura2000 can be found on related webpage of European Environment Agency. We downloaded the GIS data as shapefile in folder data/external/Natura2000_end2019_Shapefile:
path_natura2000 <- here::here(
"data",
"external",
"Natura2000_end2019_Shapefile"
)
if (!file.exists(path_natura2000)) {
temp <- tempfile()
download.file(
"https://cmshare.eea.europa.eu/s/n5L8Lrs9aYD775S/download",
temp
)
unzip(temp, exdir = path_natura2000)
unlink(temp)
}
We read the shapefile:
protected_areas <- st_read(path_natura2000)
## Reading layer `Natura2000_end2019_epsg3035' from data source `C:\Users\damiano_oldoni\Documents\INBO\repositories\pipeline\data\external\Natura2000_end2019_Shapefile' using driver `ESRI Shapefile'
## Simple feature collection with 27845 features and 6 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 746551.8 ymin: 905053.4 xmax: 6506130 ymax: 5298655
## CRS: 3035
protected_areas_crs <- st_crs(protected_areas)
Check whether grid and protected areas have same CRS:
be_grid_crs == protected_areas_crs
## [1] TRUE
We are interested in protected areas in Belgium:
protected_areas <-
protected_areas %>%
filter(MS == "BE")
Summary by site type:
protected_areas %>%
as_tibble() %>%
group_by(SITETYPE) %>%
count()
The Belgian regions are available by package BelgiumMaps.StatBel:
data(BE_ADMIN_REGION)
BE_ADMIN_REGION <- st_as_sf(BE_ADMIN_REGION)
In this section we add three columns to know whether a protected area belong to one or more Belgian regions.
Add column flanders to indicate whether the protected area is a Flemish protected areas (TRUE/FALSE):
protected_areas_flanders <- tibble(
SITECODE = c(
"BE2200043",
"BE2200032",
"BE2100017",
"BE2300006",
"BE2500121",
"BE2100024",
"BE2400008",
"BE2200030",
"BE2200037",
"BE2200041",
"BE2200042",
"BE2200038",
"BE2101437",
"BE2219312",
"BE2500831",
"BE2100323",
"BE2100045",
"BE2100020",
"BE2200028",
"BE2400009",
"BE2200033",
"BE2300044",
"BE2400012",
"BE2200626",
"BE2200727",
"BE2200029",
"BE2400010",
"BE2223316",
"BE2301235",
"BE2301336",
"BE2300007",
"BE2422315",
"BE2524317",
"BE2100026",
"BE2100040",
"BE2200031",
"BE2500003",
"BE2200034",
"BE2100424",
"BE2400014",
"BE2101538",
"BE2200039",
"BE2500932",
"BE2500001",
"BE2218311",
"BE2101639",
"BE2217310",
"BE2300222",
"BE2400011",
"BE2500002",
"BE2220313",
"BE2200035",
"BE2100015",
"BE2300005",
"BE2301134",
"BE2100016",
"BE2100019",
"BE2200036",
"BE2500004",
"BE2200525",
"BE2501033",
"BE2221314"
)) %>%
mutate(flanders = TRUE,
wallonia = FALSE,
brussels = FALSE)
Add column brussels to indicate whether the protected area is located in Brussels region (TRUE/FALSE):
protected_areas_brussels <- tibble(
SITECODE = c("BE1000001",
"BE1000002",
"BE1000003")) %>%
mutate(flanders = FALSE,
wallonia = FALSE,
brussels = TRUE)
Areas not belonging to any region as they are under federal administration:
protected_areas_federal <- tibble(
SITECODE = c("BEMNZ0001",
"BEMNZ0002",
"BEMNZ0003",
"BEMNZ0004",
"BEMNZ0005")) %>%
mutate(flanders = FALSE,
wallonia = FALSE,
brussels = FALSE)
All the other protected areas are in Wallonia (wallonia = TRUE):
protected_areas_wallonia <- tibble(
SITECODE = as.character(protected_areas$SITECODE[
!protected_areas$SITECODE %in% c(protected_areas_flanders$SITECODE,
protected_areas_brussels$SITECODE,
protected_areas_federal$SITECODE)
])) %>%
mutate(flanders = FALSE,
wallonia = TRUE,
brussels = FALSE)
We can now merge the regional information
protected_areas_regional_info <-
bind_rows(protected_areas_flanders,
protected_areas_wallonia,
protected_areas_brussels,
protected_areas_federal) %>%
mutate(SITECODE = as.factor(SITECODE))
to add it to protected_areas:
protected_areas <-
protected_areas %>%
left_join(protected_areas_regional_info,
by = "SITECODE")
Preview:
protected_areas %>%
as_tibble() %>%
select(SITECODE, SITENAME, SITETYPE, flanders, wallonia, brussels) %>%
head()
Number of protected areas per type:
protected_areas %>%
st_drop_geometry() %>%
group_by(SITETYPE) %>%
count()
In Flanders:
protected_areas %>%
filter(flanders == TRUE) %>%
st_drop_geometry() %>%
group_by(SITETYPE) %>%
count()
# Transform to wgs84
protected_areas_wgs84 <-
protected_areas %>%
st_transform(crs = 4326)
# Flemish protected areas
protected_areas_wgs84_flanders <-
protected_areas_wgs84 %>%
filter(flanders == TRUE)
prot_area_fl_palette <- colorFactor(
topo.colors(nrow(protected_areas_flanders)),
protected_areas_flanders$SITECODE)
leaflet() %>%
addTiles() %>%
addPolygons(data = protected_areas_wgs84_flanders %>%
filter(SITETYPE == "A"),
fillColor = ~prot_area_fl_palette(SITECODE),
fillOpacity = 0.7,
color = "black",
weight = 0.5,
opacity = 1,
label = ~SITENAME,
popup = ~SITECODE,
group = "type A") %>%
addPolygons(data = protected_areas_wgs84_flanders %>%
filter(SITETYPE == "B"),
fillColor = ~prot_area_fl_palette(SITECODE),
fillOpacity = 0.7,
color = "black",
weight = 0.5,
opacity = 1,
label = ~SITENAME,
popup = ~SITECODE,
group = "type B") %>%
addLayersControl(
overlayGroups = c("type A", "type B"),
options = layersControlOptions(collapsed = FALSE))
In Wallonia:
protected_areas %>%
filter(wallonia == TRUE) %>%
st_drop_geometry() %>%
group_by(SITETYPE) %>%
count()
protected_areas_wgs84_wallonia <-
protected_areas_wgs84 %>%
filter(wallonia == TRUE)
prot_area_wallonia_palette <- colorFactor(
topo.colors(nrow(protected_areas_wallonia)),
protected_areas_wallonia$SITECODE)
leaflet() %>%
addTiles() %>%
addPolygons(data = protected_areas_wgs84_wallonia %>%
filter(SITETYPE == "A"),
fillColor = ~prot_area_wallonia_palette(SITECODE),
fillOpacity = 0.7,
color = "black",
weight = 0.5,
opacity = 1,
label = ~SITENAME,
popup = ~SITECODE,
group = "type A") %>%
addPolygons(data = protected_areas_wgs84_wallonia %>%
filter(SITETYPE == "B"),
fillColor = ~prot_area_wallonia_palette(SITECODE),
fillOpacity = 0.7,
color = "black",
weight = 0.5,
opacity = 1,
label = ~SITENAME,
popup = ~SITECODE,
group = "type B") %>%
addPolygons(data = protected_areas_wgs84_wallonia %>%
filter(SITETYPE == "C"),
fillColor = ~prot_area_wallonia_palette(SITECODE),
fillOpacity = 0.7,
color = "black",
weight = 0.5,
opacity = 1,
label = ~SITENAME,
popup = ~SITECODE,
group = "type C") %>%
addLayersControl(
overlayGroups = c("type A", "type B", "type C"),
options = layersControlOptions(collapsed = FALSE))